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Abstract 

A series of weak-coupling perturbation theories which include the lowest- 
order vertex corrections are applied to the attractive Holstein model in in- 
finite dimensions. The approximations are chosen to reproduce the iterated 
perturbation theory in the limit of half-filling and large phonon frequency 
(where the Holstein model maps onto the Hubbard model). Comparison is 
made with quantum Monte Carlo solutions to test the accuracy of different 
approximation schemes. 
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I. INTRODUCTION 



The theory of superconductivity, as first developed by Bardeen, Cooper, and Schrieffer,0 
and then generahzed to strong couphng by Migdaljl Ehashberg,!, and othersill, has proven 
to be one of the most accurate theories of sohd-state physics. Properties of conventional 
low- Tc materials are typically explained to accuracies of one percent or better. Newly dis- 
covered materials, however, which have moderate Tc's, do not fit into the parameter regimes 
studied so successfully in the 50's and 60's. These materials, such as the A15 compounds, 
Bai-ajKa-BiOs, and K3C60 have large phonon energy scales relative to the inverse electronic 
density of states, so that both the effect of the energy dependence of the bare electronic den- 
sity of states (i. e., nonconstant density of states), and the effect of vertex corrections may 
become important in their description. We examine here a series of different weak-coupling 
perturbation theories that include both the effects of nonconstant density of states and of 
vertex corrections to ascertain what methods should be used for real materials calculations 
of these higher compounds. Since the most popular implementation of Migdal-Eliashberg 
theory is manifestly particle-hole symmetric (because of the neglect of the energy depen- 
dence of the electronic density of states) it is not clear what the most accurate approximation 
scheme is when the electron filling is doped away from half-filling and the energy depen- 
dence of the electronic density of states becomes important. Our strategy is to solve a model 
system (the Holstein model) in the infinite-dimensional limit via the dynamical mean field 
theory. This allows us to compare numerically exact quantum Monte Carlo solutions (in the 
thermodynamic limit) with the different perturbative approximations. 

The Holstein mo deli consists of conduction electrons that interact with local (Einstein) 
phonons: 

H = E (4cA.- + cLc,.) + T.i9X, - /i)(n,t + n^i) + ^Mfi^ Y^^'^+W^. (1) 

In Eq. (|l]), cjo- (cjo-) creates (destroys) an electron at site j with spin a, = c^j^jCja is 
the electron number operator, and Xj (pj) is the phonon coordinate (momentum) at site 
j. The hopping is chosen to be between the nearest neighbors (j and k) of a hypercubic 
lattice in d-dimensions and the unit of energy is the rescaled matrix element t* (so that 
the kinetic energy remains finite as d ^ 00). The phonon has a mass M (chosen to be 
M = 1), a frequency Q, and a spring constant k = MQ"^ associated with it. The deformation 
potential is denoted by g and it governs the strength of the coupling of electrons to phonons. 
The effective electron-electron interaction strength (mediated by the phonons) is then the 
bipolaron binding energy 

A chemical potential /i is employed to adjust the total electron filling with n = U corre- 
sponding to half-filling in the exact solution. 

In the instantaneous limit where U remains finite and g and VL are large compared to the 
bandwidth {g,VL —>■ 00, U = finite), the Holstein model maps onto the attractive Hubbard 
modeli 
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with f/ defined by Eq. (0). 

Two cases of the Holstein Hamiltonian have well-estabhshed perturbative expansions. 
In the hmit where the phonon frequency becomes small ^ 0, but U remains finite, 
Migdal-Eliashberg theorj^'i is known to be an accurate approximation. Migdal-Eliashberg 
theory is a self-consistent Hartree-Fock approximation that employs fully dressed phonon 
propagators and neglects all vertex corrections. Typically, the energy dependence of the 
electronic density of states is also neglected, so the theory is evaluated with a constant 
density of states. (This latter approximation is always particle-hole symmetric and maps 
onto the limit of half-filling when the nonconstant density of states is used.) The second 
limit is the large phonon frequency limit Q oo with U also remaining finite. In this case, 
a truncated perturbative expansion through second orderi'i (which includes the lowest-order 
vertex correction) is known to be accurate at half-filling for a large range of U values, because 
it properly reproduces both the weak-coupling and strong-coupling limits of the Hubbard 
model. It would be nice to construct an approximation scheme that continuously connects 
these two limits as the phonon frequency is varied. However, no simple approximation can 
be found, because the set of diagrams that dress the phonon propagator in the small-f2 limit 
do not correctly dress the phonon propagator in the large- limit; in the large-frequency 
limit the interaction is only between up and down spins, but in the small-frequency limit the 
interaction is between all spins. But, if one is willing to examine perturbative expansions 
that are truncated with respect to the fiuctuating diagrams, then approximations can be 
constructed that agree with the two known limits through the order of diagrams included 
in the expansion. The description of these different methods is both subtle and technical, 
and will be covered in detail in Section II. Here we only want to comment that the previous 
calculations for the electron-phonon problem that were called the iterated perturbation 
theory! (IPT) are actually based on a truncated perturbation theory about the Hartree 
mean-field theory solution, which does not reproduce the Migdal-Eliashberg (Hartree-Fock) 
limit properly. A more promising approach for all Q is to construct a truncated perturbation 
theory about the Hartree-Fock mean-field theory solution, which is done here. In addition, 
we also examine some simple methods that can be used to repair inconsistencies that develop 
in the IPT as the system is doped off of half-filling. This method entails self-consistently 
renormalizing the higher-order fiuctuations of the Hartree or Fock diagrams to all orders, 
which allows for the electron filling on the lattice to be different from the electron filling of the 
mean-field theory solution when the higher-order fiuctuations are included. Unfortunately, 
we find that all attempts to produce an accurate perturbation theory that reproduces the 
IPT in the instantaneous limit are not very accurate for moderate phonon frequencies. 

These perturbative approaches, and the many-body problem in general, simplify when 
the infinite-dimensional limitty is taken {d oo). Then the many-body problem becomes a 
local problem that retains its complicated dynamics in time. The hopping integral is scaled 
to zero in such a fashion that the free-electron kinetic energy remains finite while the self 
energy for the single-particle Green's function and the irreducible vertex functions have no 
momentum dependence and are functionals of the local Green's function0^ii. This limit 
retains the strong-correlation effects that arise from trying to simultaneously minimize both 
the kinetic energy and the potential energy, and hence has relevance for three-dimensional 
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materials. 

Of course, we can also solve the infinite-dimensional Holstein model using a quantum 
Monte Carlo method, which contains all effects due to phonon renormalizations, vertex cor- 
rections, and nonconstant electronic density of states. We use these solutions as a benchmark 
to test the accuracy of the different approximation methods and to determine what is the 
most fruitful approximation for the electron-phonon problem. 

We employ a Green's function formalism to solve the many-body problem. The local 
Green's function is defined to be 

GiocituJn) = - dre t^^^^W) ' ^^'> 

and is calculated directly from a self-consistent quantum Monte Carlo procedure described 
elsewhere.lli Static two-particle properties can also be determined since the irreducible vertex 
function is locaH The static susceptibility for charge- density- wave (CDW) order is given 
by 



CDW 



;q) = E e^'^-(^^-^^)T f'dr f' dr'[{n,.{T)n,Ar')) - (n,.(r))K.,(r'))], 

2^ R,-R,aa' ^0 Jo 

= rEx^r(q), (5) 
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at each ordering wavevector q [the indices m and n denote Matsubara frequencies iun = 
i7vT{2n + 1)]. Dyson's equation for the two-particle Green's function becomedllE 

x^r (q) = xli^)5mn - rEx^(q)r^rC'''^(q)' (6) 

p 

with r^^'^ the irreducible vertex function in the CDW channel. 

The bare CDW susceptibility Xn(q) iii Eq. (H) is defined in terms of the dressed single- 
particle Green's function 

Xn(q) = -^EG^n(k)Gn(k + q) 
k 

2 2 

1 r°° e^" 

TT 

with all of the wavevector dependence described by the scalarlll' 

^ ^(Ol = E?=iCOsq,/c?. 

The integral for ^^X) in Eq. (J^) can then be performed analyticallyt^ for the "checker- 
board" CDW phase Xn(-^ = ~1) = ~Gnlii^n + A* ~ The irreducible vertex function 
^mn^ is either directly calculated in a perturbative expansion (described below) or is deter- 
mined by inverting the Dyson equation in Eq. (|^) for the local charge susceptibility (QMC). 

A similar procedure is used for the singlet s-wave superconducting (SC) channel. The 
corresponding definitions are as follows: The static susceptibility in the superconducting 
channel is defined to be 



dy-. — / dz . = (7) 

oo iu;n + f^-2^n-y J-c^ iuJn + l^-^n- X{q)y - zJl- X^{q) 



q) = ^ E ^^''-^'''-'''^T f' dr f' dr\c,,{T)c,,{rH^^^^ (8) 

Rj-Rfe ° ° mn 
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for superconducting pairs that carry momentum q; Dyson's equation becomes 



P 

with r^*^ the corresponding irreducible vertex function to describe the SC channel; the bare 
pair-field susceptibility becomes 

Xn'(q) = ^EGn(k)G'-n-l(-k + q) 
k 

2 2 

1 poo p~y roQ p—^ 

= - / dy / dz , (10) 

ttJ-oo tuJn + fi-J:n-y J-oo icj_„_i + /i - - X(q)y - 2^1 - X2(q) 

with the special value Xn'i-^ — ^) — ~lHiG'„/Im(ici;„ — S„) for the SC pair that carries no net 
momentum; and finally the irreducible vertex function is also either directly calculated in a 
perturbative expansion (described below) or is determined by inverting the Dyson equation 
in Eq. (|^) for the local pair-field susceptibility (QMC). 

The transition temperature of the infinite-dimensional Holstein model is then found by 
calculating the temperature at which the relevant susceptibility diverges (CDW or SC). This 
transition temperature is found by locating the temperature where the scattering matrix (in 
the relevant channel) 

has unit eigenvalueS (note that the local Green's functions are always used in the evaluation 
of the bare susceptibility xo)- 

The remainder of the paper is arranged as follows: In Section II we describe the formalism 
employed to derive the different truncated perturbative expansions for both the electronic 
self energy and for the irreducible vertex functions. Particular attention is paid to the 
mapping from the lattice problem onto an effective impurity problem, and how to extract 
a perturbative expansion that includes the lowest-order fluctuations beyond a mean-field 
theory solution. In Section III we describe the details of the calculational procedure and 
present our results in Section IV. Section V contains a summary and our conclusions. 



II. LATTICE-IMPURITY MAPPING 

The properties of the Holstein model are calculated by mapping the lattice problem 
onto the single-impurity Wolff-Holstein mo delllii, and solving the impurity self energy by a 
truncated perturbation expansion. To define the mapping, we begin with the self-energy of 
the Holstein model, which is momentum independent in infinite dimensions, and write the 
local Green's function for the lattice as, 

Giociu;) = Y: : \ y (12) 

Y ^ - (ep - /^) - S(a;) 

where ep is the noninteracting electronic bandstructure for the lattice (assuming spin de- 
generacy = = n/2). The local self-energy S(a;) is defined by the functional, 



5 



S(c<j) = S[G'ioc]7 where S [G/oc] represents the sum of all the skeleton diagrams gener- 
ated by the perturbation theory, with U as the expansion parameter (in a Baym-Kadanoff 
expansionlll) . On the other hand, the same set of skeleton diagrams appears in the self- 
energy of an impurity problem described by the Wolff-Holstein Hamiltonian, 

Himp = Ho + gxoino^ + hq^) + ^ + ^Mil^xl (13) 

where Hq describes a band of non-interacting electrons, Hq^ is the number operator for 
conduction electrons of spin a at the impurity site and xq (po) is the "impurity" phonon 
coordinate (momentum). The renormalized electron propagator at the impurity site can be 
written as, 

G.„p(^) = , ^ ■ (14) 

Go [uj) - S(u;) 

Here Gq is the free-electron propagator at the impurity site (which is often called the 
effective- medium propagator), and S(a;) describes the effects due to the coupling to the 
local phonon. Note that the effective-medium Gq is not equal to the noninteracting lattice 
propagator except when U = 0. 

Since the impurity and the lattice self-energy functionals, E [Gj^p] and S [Gioc] are the 
same, the self energies will also be the same if Gimp{uj) = Giodyj). Thus, the lattice problem 
maps onto the impurity problem, provided the effective medium propagator Gq is adjusted so 
that the right hand sides of Eqs. ([T2|) and (14) are equal. In general, the self-energy functional 



S [Gioc\ is not known. But equivalent expansions can also be made for the impurity self- 
energy by rearranging the skeleton-diagram expansion. For example, a functional defined on 
the effective-medium propagator Gq{uj) can be employed, such that S(ci;) = Sq [Gq], where 
Eq [Go] represents the sum of all connected graphs generated by Wick's theorem (without 
any resummations) . If the graphs for S [G^oc] and Sq [Gq] are summed to all orders, then 
the self energies must agree. Furthermore, the effective medium Gq{uj) can be uniquely 
determined by either Gq\uj) = GrolH + S [G] or Gq^uj) = G,~J(cj) + Sq [Gq]. 

Since the exact self-energy for the Wolff-Holstein model with an arbitrary density of states 
is not known, our strategy is to make a truncated expansion for the impurity self-energy 
functional S^pp [Gq], and define the self-energy for the lattice problem via S(ci;) = S^pp [Gq], 
The self consistency is then imposed on Go(cc;), so that Giod^^) = Gimp{oj). The hope is, 
that a controlled expansion of the impurity self-energy will lead to an accurate solution of 
the impurity problem, and of the resultant lattice problem. Here, we use the lowest-order 
Yosida-Yamada expansion, which is known to provide reliable answers for the Wolff model 
describing magnetic impurities00. 

The Yosida-Yamada expansion for the effective impurity problem (with an unknown den- 
sity of states) is obtained by a partial resummation of diagrams that allows us to construct 
the self-energy functional in terms of diagrams that involve either Hartree or Hartree-Fock 
Green's functions. That is, one considers the self-energy corrections with respect to the 
mean-field solution, and rewrites the Dyson equation as, 

Gimpiuj) = J— 7~Tr:i -pMFir^ ]■ i^^) 
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where Eyy is the sum of all self-energy diagrams defined in terms of the mean-field propa- 
gator for the impurity problem, Gmf, and 




(16) 



It is clear that functionals S^y [Gmf] defined with different mean-field propagators are not 
the same, and it is not known a priori which truncated approximation for the impurity self 
energy leads to the most accurate result for the lattice problem. The mean-field solution 
has an average electron filling umf defined by 



with f{w) = 1/[1 -|- exp{i3uj)] the Fermi factor (the factor of 2 is from the spin summation). 
This filling is usually different from the electron filling on the lattice which is determined 
by the same equation, but with Gmf replaced by Gioc- In this contribution, we consider 
four different expansions to the self energy and compare the corresponding results to the 
quantum Monte Carlo solution. 

Our first approximation is called the Hartree expansion (H) and it includes a truncated 
perturbation-theory expansion about the Hartree mean-field solution, including all diagrams 
through second order. This expansion should not be confused with the Hartree approxima- 
tion, which does not include any of the second-order vertex corrections. These self-energy 
corrections with respect to the Hartree mean-field solution are constructed by using the 
Hartree self-energy diagram shown in Figure 1, to construct the mean- field propagator of 
Eq. (pISl). That is, we take 



which implies that the Hartree self-energy insertions are included to all orders in the pertur- 
bative expansion. Here is the (Hartree) mean- field particle number defined by Eq. (p!7|) 
with Gmf = Gh- Next, the Hartree mean-field Green's function is used in evaluating the 
truncated Yosida-Yamada self-energy functional to determine the approximate self-energy. 
This functional is given in Figure 1 through second-order in U where the electronic propa- 
gators (solid lines) are the Hartree mean-field propagators from Eq. ([I6|), the phonon propa- 
gators (wiggly lines) are bare D{uj) = — l/[M(fi^ — cu^)], and the impurity self-energy which 
renormalizes the effective field Gq in Eqn.([l^) becomes 



Note that this procedure is summing an infinite class of diagrams (the Hartree self-energy 
insertions), but is truncating the perturbation theory with respect to the Hartree mean- 
field-theory solution to include both the Fock and the second-order fiuctuating terms. It 



approximation^, since the renormalization of the effective medium by the Hartree self energy 
(to construct the Hartree mean-field Green's function) can be absorbed by a redefinition of 
the chemical potential (since the Hartree self energy just provides a frequency-independent 
shift), with the exception being the inclusion of the fifth diagram, which renormalizes the 
Hartree diagram by the Fock self-energy insertion. The omission of the fifth diagram in 




(17) 



Smf — — Uuh, 



(18) 



j:{uj) = j:h + ^yy[Gh]. 



(19) 
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the older IPT papem may be thought to be innocuous, because that diagram vanishes 
at half-filhng and should not affect the results much off of half filling. However, when 
the coupling strength becomes large, it's effects do become strong, as shown in the next 
section. To reiterate, the difference between the present expansion and the older IPT work 
is that the fluctuating diagrams in Syy[G//] are evaluated with Gh which does not include 
the frequency-independent shift from the fifth diagram of Figure 1. The IPT calculation 
included the frequency-independent shifts to all orders, and hence evaluated the first four 
diagrams of Figure 1 using a different mean- field propagator than Gh- The inclusion of 
the Fock self-energy insertion into the Hartree diagram is just one of the subtle, and often 
neglected diagrams that needs to be included in a truncated approximation that includes 
all diagrams up to a given order. 

Our second approximation is called the Hartree-Fock expansion (HF) and is obtained by 
using the Hartree-Fock self-energy 

Smf = ^HF = Uuhf + (20) 

to define the mean-field propagator in Eqn.(p!6D, and including all second-order diagrams 
with respect to the Hartree-Fock mean-field solution (once again, this should not be confused 
with the Hartree-Fock approximation which does not include any of the second-order vertex 
corrections). The (Hartree-Fock) mean-field particle number uhf is calculated by using 
G/fi? in Eqn.(0), and the Fock self energy J^p satisfies 

j:piuj) = g'l deGHFiuj - t)D{t)f{t), (21) 

with -D(e) the bare phonon propagator again. This mean- field solution sums both the Hartree 
and Fock self-energy insertions to all orders, and would be identical to Migdal-Eliashberg 
theory if the dressed phonon propagator was employed in Eq. (|2T|) rather than the bare 
propagator. As described in the introduction, we are forced to use the bare propagator if we 
want to reproduce the IPT limit in the large-phonon-frequency limit. The Yosida-Yamada 
functional defined on Ghf is given by the diagrams in Figure 2 where the solid lines are 
now Hartree-Fock propagators. The self-energy correction to Gq for the HF expansion can 
then be written as 

1:{uj) = 1:hf + ^^^[Ghf\. (22) 

The self-energy diagrams given in Figs. 1 and 2 are evaluated on the imaginary axis 
by using the standard rules {explicit expressions are given in Eq. (16) of Ref. |^ [and the 
suitable modifications for the Hartree-Fock case]}, and the Matsubara summations occurring 
in T^yyI'^h] and Eyy [Gj^f] are performed numerically for both the Hartree and the Hartree- 
Fock expansions. 

If the perturbative expansion for the self-energy is extended to higher-order in U we find 
additional frequency-independent diagrams [which sum up to U{n — nMF)] in addition to the 
frequency-dependent diagrams (similar to the fifth diagram in Figure 1). Here n is the fully 
renormalized particle number of the lattice which is calculated by Eq. ( p!7D with Gioc replacing 
Gmf- These diagrams arise simply from the fact that in the exact skeleton expansion the 
Hartree diagram is evaluated with the fully dressed Green's function (yielding Un) rather 
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than with the mean- field Green's function (which yields Uumf)- At half- filling, the two 
fillings n and umf are usually equal, but they need not be the same (and never are) away 
from half filling. We introduce a new self-energy functional Syy[G//] (that is based on the 
Yosida-Yamada functionals defined above) that has all of the higher-order renormalizations 
of the Hartree diagram removed from it [as shown, through third-order in Figure 3 (a)]. 
Then we have an exact relation S(co') = Un + J^yyIGh] for the self energy S(a;) defined in 
Eq. ([T^). In this contribution, the two self-energy functionals Syy and Syy differ by one 
diagram (the fifth diagram in Figure 1). We can employ this exact relation to formulate 
the renormalized-Hartree expansion^ (also called the n-consistent approximation), which 
enforces this self-consistency condition on the renormalized particle number. That is, for a 
given /i, g, uh and Syy(ci;), we solve the transcendental equation, 



n 



--ImY^I doo ^- , (23) 



to determine the impurity electron filling n, and to obtain the correction to the real part 
of the total self energy. For the Anderson and Wolff impurity models, and for the two- 
dimensional Hubbard model, this renormalized-Hartree approach (RH) expands the re- 
gion of the validity of the perturbation theory, and allows for the description of stronger 
correlationsilii. The hope is that this renormalized-Hartree scheme will repair some of 
the inaccuracies of the truncated perturbation theories away from half filling. It is a much 
simpler approximation to study than a scheme that tries to interpolate between the weak- 
coupling and atomic limitH (which is significantly more challenging to formulate for the 
electron-phonon problem) . 

A similar approach can be made for summing higher-order diagrams in the Hartree- 
Fock expansion. In this case, however, an n-consistent approximation, that sums only the 
frequency-independent diagrams to all orders, is likely to be less accurate than one that 
sums the renormalizations based on both the Hartree and Fock diagrams. Hence, we form 
the exact relation for the self energy of Eq. (|T3) 



Un + g^l deGUuj - e)D{e)f{e) + (24) 



and employ it to evaluate the renormalized-Hartree-Fock expansion (RHF) where the func- 
tional T,yy{Ghf) is truncated at second order [see Figure 3(b)]. This approximation is 
formed by evaluating both the Hartree and Fock diagrams with the local Green's function 
(instead of Ghf) but the second-order diagrams continue to be evaluated with Ghf- In 
this sense it is "halfway" between a conserving approximation (which evaluates all diagrams 
with Gioc) and the truncated expansions described above. 

The only remaining objects that need to be determined are the irreducible vertices for 
the CDW and SC instabilities. These vertices are calculated with the mean-field Green's 
function (either Hartree or Hartree-Fock) and have an identical functional form for all dif- 
ferent approximation schemes. The diagrams included are presented in Figure 4(a) for the 
CDW vertex and Figure 4(b) for the SC vertex. The solid lines are Gmf {Gh or Ghf) and 
the wiggly lines are the bare phonon propagator. Explicit expressions for these diagrams on 
the imaginary axis are given in Eqs. (17) and (18) of Ref. |] with Go replaced by Gmf- 
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III. CALCULATIONAL METHODS 



The calculational methods used are straightforward. The perturbation expansion is 
carried out on the imaginary axis. We employ an iterational algorithm to solve the self- 
consistent equations, which is summarized in Figure 5: 

1. Start with an initial self-energy T,{uj); 

2. Use Eqn.(|l2D to calculate Gioc', 

3. Calculate the effective-medium propagator Gq using Eqn.(|T^ with Gimp = Gioc', 

4. Calculate the mean-field propagator using Eqn.(|l6|); 

5. Calculate the Yosida-Yamada self-energy diagrams as defined by Fig. 1 or Fig. 2, with 
and without the renormalization of the first-order diagrams (as shown in Fig. 3), to 
determine the total self-energy; 

6. If the calculation for the self-energy has not converged, then adjust the chemical po- 
tential to produce the target electron filling; 

7. Repeat steps 2-6 until the calculation has converged; 

8. Once converged, calculate the irreducible vertex functions from Fig. 4 and determine 
the maximal eigenvalue of the scattering matrix in both the CDW and SC channels 
from Eq. (|Tl]); 

9. Repeat 1-8 for another temperature T until the transition temperature Tc is found. 

We use an energy cutoff of 256 Matsubara frequencies, which provides accurate results for 
temperatures T larger than O.Olt*, and our convergence criterion is that the Green's functions 
do not change by more than one part in 10^ from one iteration to the next. We find that 
the perturbation theory typically converges in approximately 100 iterations of the main ring 
in Figure 5. Sometimes the iterated equations develop limit cycles, whose oscillations can 
be suppressed by employing standard damping methods that average the i — 1st and ith 
iterations to produce the starting point for the next iteration. We also found that if the 
coupling strength is large enough, then some of the approximate theories will have multiple 
solutions for n{fi) near half filling. In this case, the symmetry point, with fi = U becomes 
unstable (i. e., the compressibility is negative), and the Green's functions are no longer 
purely imaginary when the electron filling is half-filled. This latter result is an indication of 
the breakdown of the approximation method at such a large coupling strength. 

There also are sometechnical details that need to be discussed about the quantum 
Monte Carlo simulationsta. The results for the transition temperatures were calculated with 
a At = 0.4, and sometimes also with a At = 0.2 and then extrapolated to At = 0, when the 
correlations were large enough that the Trotter error was noticeable. The self-energies and 
irreducible vertex functions were calculated with a fixed At = 0.1 to ensure a high accuracy. 
In general, there are three sources of error to the quantum Monte Carlo simulation: (i) 
statistical error; (ii) iteration error; and (iii) systematic error. The statistical error is the 
easiest error to control, and is the smallest of the three sources of error. The iteration error 



10 



arises from performing calculations with the wrong dynamical mean-field Go, because the 
calculation has not yet converged fully. We typically iterate 8 or 9 times, which provides 
good convergence for the iterations. The systematic error is more difficult to control. It 
arises from the Trotter error, and from other (potentially unknown) sources. One surprising 
source of error arose from the choice of global updating moves for the phonon coordinate. 
Global moves that shift the phonon coordinate at every time slice by an amount randomly 
chosen between —dx/2 and dx/2 were supplemented by global shifts by an amount ranging 
between -^2g/{mVL^) — dx/2 and and -^2g / {mVt^) + dx/2. These latter moves were chosen 
to allow the phonon coordinate to shift between the two minima of the effective phonon 
potential (separated by 2g /[mVt^])., which form when the system enters the strong-coupling 
region, were preformed pairs exist above T^. When one is in the weak-coupling regime, where 
the effective phonon potential possesses a single-well structure, we find that the self-energy 
can differ by a few percent depending on whether the global moves are all chosen randomly, 
or if the supplemental global shifts are chosen to move the phonon coordinate to the "other" 
well (i.e. if the phonon coordinates lie in the right well then the global shift is chosen to 
move the system to the left, and vice versa). The accuracy was best when the global moves 
were chosen completely randomly, which is the method used here. These results differ by 
a few percent from those shown previouslyl^l, where the global moves were coupled to the 
values of the phonon coordinates, and the Trotter time slice was larger. 



IV. RESULTS AND DISCUSSION 

The parameters of the Hamiltonian in Eq. (|^) are chosen to agree with previous theoret- 
ical work, and to represent a parameter range where the vertex corrections can cause large 
effects. As such, we choose the phonon frequency to be on the order of one tenth of the band- 
width VL = 0.5t* (with the effective bandwidth of the Gaussian being about two standard 
deviations above and below the center, or approximately At*). Most of our calculations are 
performed for three different values of the interaction strength (i) g = 0.4 {U = —0.64, one 
sixth of the bandwidth) which corresponds to a fairly weak interaction where there are no 
preformed pairs, but the transition temperatures are large enough that the phase diagram 
can be determined reliably; (ii) g = 0.5 {U = —1.0, one fourth of the bandwidth) which 
is in the moderate interaction regime; and (iii) g = 0.625 {U = —1.5625, two-fifths of the 
bandwidth) which is where the system enters the strongly coupled regime, where numerous 
preformed pairs are present above Tg. 

The quasiparticle renormalization factor, the irreducible CDW vertex function, and the 
irreducible SC vertex function are plotted in Figs. 6(a-c). The vertex functions are averaged 
to show just the symmetric frequency component, since that is all that contributes to the 
eigenvalue of the scattering matrix for the CDW at half filling, or for the SC at any filling 
(this is so is because the eigenvector with maximal eigenvalue is symmetric with respect 
to a change in sign of the Matsubara frequency). The coupling strength is g = 0.4 and 
the filling is n = 1.0, which yields a bare electron-phonon coupling strength of = 0.64, 
or A = p(0)|f/| = 0.36, which would naively be viewed as quite weak coupling. However, 
this is not the case, as one can estimate the renormalized value of A from Z{0) — 1 which 
lies at about 0.7. These calculations indicate that higher-order diagrams, such as those 
that fully dress the phonon propagator, are important even when the bare coupling has 
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A = 0.36. What is interesting, is that the perturbative results seem to have the right shape, 
but need to have the frequency axis moved to the right, to hne up with the Monte Carlo 
data. Furthermore, since the perturbation theory underestimates both the renormalization 
factor (which will enhance T^) and the vertex function (which will reduce T^) one might 
expect these results to tend to cancel each other out when calculating a phase diagram, 
which is indeed what we find below. Figs. 7(a-c) show the same results for the stronger 
coupling case g = 0.5 and n = 1.0 {\U\ = 1.0 and A = 0.56). We can see that the error 
in the renormalization factor has grown, even though all four perturbative methods yield 
similar values. The CDW vertex becomes more attractive at moderate coupling, which is 
missed by the perturbation theory. The SC vertex has the right qualitative shape, but is 
underestimated at low frequencies. When the coupling strength is increased to g = 0.625 
(not shown) the situation becomes even worse, with the same qualitative features seen, a 
growing difference in the renormalization factors, and an underestimation of the size of the 
vertices. 

We also include, in Fig. 8, a result at quarter filling n = 0.5 and g = 0.4 at the same 
temperature T = 1/16 as in Fig. 6. One might have thought that the quarter-filled case 
would be approximated much better by the perturbation theory because the effective cou- 
pling strength is lower due to the reduction of the electronic density of states at the Fermi 
level, but we see that the shape of the curves is both qualitatively and quantitatively similar 
for Figs. 6 and 8. The quantitative agreement has improved, but the improvement is not 
dramatic. What is surprising is that the calculated transition temperature agrees much bet- 
ter with the QMC here than at half filling, because of the cancellation of the errors to both 
the self energy and the vertex function. Note that we only plot the SC vertex here, because 
an asymmetric average is needed for the physically relevant piece of the CDW vertex. 

In addition to calculating the single-particle and two-particle self energies, we can also 
investigate the phase diagrams of the Holstein model. Based on our results for the self 
energy and the vertex functions, it does not appear likely that any approximation will be too 
accurate for the phase diagrams, but as theorized above, it is possible that the inaccuracies 
will cancel out to produce more accurate results for the transition temperatures. Such 
behavior should be viewed as "accidental" agreement with the exact result. 

We begin by calculating the transition temperature to the commensurate (checkerboard! 
CDW insulator at half-filling as a function of the interaction strength. The QMC solutionli^ 
showed that the maximal Tc was on the order of l/25th of the bandwidth, occurring when U 
was about two- fifths of the bandwidth. The results for each of the four expansions shown here 
is presented in Figure 9. The horizontal axis is the interaction strength and the vertical axis 
is the transition temperature. One can see that the peak in versus U is not produced by 
the perturbation theory, rather the transition temperature continues to increase. This differs 
from what happened in the so-called IPT approximation where the approximate curve for Tc 
did show a peak.§ The Hartree and renormalized-Hartree expansions agree with each other 
and the IPT approximation until U is approximately —1, where they differ from each other 
because the thermodynamically stable H and RH expansions are no longer at the particle- 
hole symmetric point with fi = U, but rather have different values of n (and no longer have 
purely imaginary Green's functions when evaluated along the imaginary axis). This is the 
point where the kink appears in the Tc{g) curve. The remarkable agreement of the IPT 
approximation in predicting both the peak position and the peak height accurately is really 
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just a coincidence, since neither the self energy nor the vertex function is reproduced well 
at that value of the interaction strength (see Figures 6-8). There is a very delicate balance 
between the self-energy terms (which tend to reduce as the interaction increases) and the 
vertex-function terms (which tend to increase as the interaction increases) that causes 
to have a peak for the IPT approximation, but continue to increase for the approximations 
in this contribution. 

We end with a calculation of the phase diagrams off of half-filling, which show the tran- 
sition from the checkerboard CDW insulator to a SC. These results are plotted in Figure 10. 
The horizontal axis is the electron filling and the vertical axis is the transition temperature 
to either a CDW or SC phase. The CDW-SC phase boundary lies at the points where the 
phase diagrams have a slope discontinuity. Three values of interaction strength are included, 
and both the approximate solutions and the QMC results are presented (the lines through 
the QMC data are just a guide to the eye). The agreement for is quite good for the 
weak-coupling value g = 0.4, with all approximations going through the QMC data for the 
SC phase, and overestimating only slightly for the CDW phase. Here, the Hartree-Fock 
expansion is actually worse than the Hartree expansion for the CDW phase. The calculation 
also seems to underestimate the location of the CDW-SC phase boundary. The case with 
g = 0.5 begins to show how the solutions develop qualitatively wrong behavior near half 
filling. All but the Hartree-Fock expansion show a suppression of as the filling approaches 
1. Such a suppression was never seen in the QMC simulations, and is likely an artifact of 
the non-conserving nature ot the approximation J. It continues to be clear that the SC solu- 
tions are approximated better than the CDW solutions, and that the approximations behave 
worst near half-filling. In particular, it is important to note that the simple n-consistent 
scheme, which worked so well for the Hubbard model, does not appear to work as well for 
the electron-phonon problem examined here, because the phase diagrams (H and RH) con- 
tinue to have qualitatively incorrect behavior. Finally we examine the strongly coupled case 
g = 0.625. Here the renormalized expansions do a better job at repairing the qualitatively 
wrong behavior of a Tc suppression near half-filling, but they do so by greatly enhancing the 
Tc's at and near half-filling from the truncated expansions. It is clear, that at this large a 
value of the coupling strength, none of these weak-coupling approaches is doing a good job 
at predicting the phase diagrams. It is interesting to note that the close agreement of the 
Hartree and renormalized Hartree calculations shows that the majority of the shift of the 
chemical potential comes from the Fock dressing of the Hartree diagram. Surprisingly, the 
inclusion of this diagram has a large effect on the phase diagram when \U\ > 1. 

V. CONCLUSIONS 

What are the conclusions that can be drawn from this work? Even at a weak value of the 
bare interaction strength {g = 0.4) both the self energy and the irreducible vertex functions 
are not approximated accurately by the perturbation theory. This implies that third-order 
(and higher-order) diagrams are important and cannot be neglected. For example, it is 
possible that a scheme that uses an expansion with dressed phonon propagators would 
be more accurate, which is work in progress. Nevertheless, the perturbation theory does 
appear to produce a more accurate approximation to the phase diagrams than the individual 
self-energies. On the other hand, these results show that there is no simple scheme that 
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will allow one to approximate the electron-phonon problem accurately for all values of the 
phonon frequency. So one should not necessarily rely on approximations that employ purely 
electronic models (such as the Hubbard model) to carry over to the electron-phonon problem 
and work well when the phonon frequency is small to moderate. Instead, it is more fruitful 
to work on generalizations of the Migdal-Ehashberg theory that work with dressed phonons, 
but include higher-order nonadiabatic effects such as vertex corrections. This should save a 
lot of time in formulating accurate approximation methods for real materials that have large 
enough phonon energy scales that the effects of nonconstant density of states and of vertex 
corrections cannot be neglected. Working with dressed phonons actually makes a lot of sense 
for the electron-phonon problem, because experiments can directly measure the dressed 
phonon spectral function, so it is readily available for use in real materials calculations. 
Furthermore, the dressed phonon propagator can be extracted from the QMC simulation 
and used in the perturbation theory for the electrons, just as is done from experiments on 
real materials. Work along these lines is currently in progress. 
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FIGURES 



FIG. 1. Feynman diagrams included in the Hartree expansion. The top figure shows the Hartree 
diagram, which is incorporated into Gmf, and the bottom shows the Yosida-Yamada functional 
about the Hartree mean-field solution expanded through second order. The wiggly lines are bare 
phonon propagators, and the straight lines are Hartree (mean-field) propagators. The fifth diagram, 
that involves the Fock dressing of the Hartree diagram was not included in the previous work. 

FIG. 2. Feynman diagrams included in the Hartree-Fock expansion. The top figure shows 
the Hartree and Fock diagrams, which are incorporated into Gmf, and the bottom shows the 

Yosida-Yamada functional about the Hartree-Fock mean-field solution expanded through second 
order. The wiggly lines are bare phonon propagators and the straight lines are Hartree-Fock 
(mean-field) propagators. 

FIG. 3. Modified Yosida-Yamada functional for the renormalized expansions, (a) The renor- 
malized-Hartree expansion (also known as the n-consistent approximation), in which all dressings 
of the Hartree diagram are removed from the Yosida-Yamada functional, (b) The renormalized 

Hartree-Fock expansion, in which all dressings of the Hartree and the Fock diagrams arc removed 
from the Yosida-Yamada functional. The wiggly lines are the phonon propagators and the straight 
lines are the corresponding mean-field Green's functions (Hartree or Hartree-Fock). 

FIG. 4. Irreducible vertex functions for (a) charge-density-wave and (b) superconducting order. 
The wiggly lines are the bare phonon propagators and the straight lines are the corresponding 
mean-field propagators (Hartree or Hartree-Fock). 

FIG. 5. Iterative algorithm for solving the self-consistent perturbation theory. As described 
in the text, one starts from an initial self-energy, then constructs the local Green's function, the 
effective-medium Green's function, and then the mean-field Green's function. The self-energy is 
then determined from the Yosida-Yamada expansion, and if the calculation is not converged, then 
the chemical potential is updated to adjust the electron-filling and the process is repeated. If the 
calculation has converged, then the irreducible vertex functions are determined from Figure 4, and 
the largest eigenvalue of the scattering matrix is then calculated. 

FIG. 6. Self energy and irreducible vertex functions for the Holstein model at half filling 
(n = 1.0) and weak-coupling g = 0.4. The temperature is just above Tc at T = 1/16. Fig- 
ure 6(a) shows the quasiparticle renormalization factor minus 1 as a function of frequency. The 
four different expansion schemes (Hartree, Renormalized Hartree, Hartree-Fock, and Renormalized 
Hartree-Fock) arc plotted with lines, and the quantum Monte Carlo is plotted with solid dots. We 
believe the combined errors on the simulation to be on the order of a few percent. Figure 6(b) is 
the CDW vertex and Figure 6(c) is the SC vertex. 
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FIG. 7. Self energy and irreducible vertex functions for the Holstein model at half filling 
{n = 1.0) and moderate-coupling g = 0.5. The temperature is just above at T = 1/9. Figures 
7(a), (b), and (c) are the quasiparticle renormalization factor minus 1, the CDW vertex, and the SC 
vertex respectively. Note how the agreement with the perturbation theory worsens as the coupling 
strength is increased. 

FIG. 8. Self energy and irreducible vertex functions for the Holstein model at quarter filling 
{n = 0.5) and weak-coupling g = 0.4. The temperature is the same as in Figure 6, T = 1/16. 
Figure 8(a) is the quasiparticle renormalization factor minus 1 and Figure 8(b) is the SC vertex. 
Note how the perturbation theory has improved, but less than what would have been expected 
naively. 

FIG. 9. Phase diagram of the Holstein model at half filling. The four different expansion 
schemes (lines) are compared to the quantum Monte Carlo (dots). Note how all perturbative 
approximations do not show a peak in the CDW transition temperature at half filling, but rather 
continue to increase. The Hartree and Renormalized Hartree expansions are identical until the 
solution at half filling (with chemical potential n = U) becomes unstable for the RH expansion, 
and the Green's functions acquire real parts. 

FIG. 10. Phase diagram of the Holstein model away from half filling for three different coupling 

strengths g = 0.4, g = 0.5, and g = 0.625. The lines through the QMC data points are a guide 
to the eye. Note that the SC transition temperatures arc approximated better than the CDW 
transition temperatures, and that the phase diagrams are predicted more accurately than would 
be expected from the individual self-energies or vertex functions. As the coupling strength is 
increased the T^s are all enhanced significantly, and the location of the CDW-SC phase boundary 
is predicted less accurately. The renormalized expansions do sometimes improve the qualitative 
shape of the phase boundaries, to show a maximum Tc at half filling, but the perturbative approach 
is definitely breaking down as the coupling strength becomes larger than 0.5. 
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Figure 1 , Phys. Rev. B, Freericks, et al. 




Figure 2, Pliys. Rev. B, Freericks, et al. 
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Figure 3, Phys. Rev. B, Freericks, et al. 
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Figure 4, Phys. Rev. B, Freericks, et al. 
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Figure 5, Phys. Rev. B, Freericks, et al. 



21 




Freericks, et al., Phys. Rev. B, Fig. 6 (a) 



22 







-0.5 - 



-1 - 



1 


1 1 1 

(b) 




II 

H 




n> 1 1 

RH 




!_! U 

Hr 




RHF 




• MC" 






1 


• • • • 

1 1 1 







12 3 4 5 

Matsubara Frequency ico^ 

Freericks, et aL, Phys. Rev. B, Fig. 6 (b) 



23 



C\l 







- -0.2 - 



o 

CO 



-0.4 I- 
+ -0.6 
^. -0.8 h 



o 
if) 



-1 







1 


(c) 




H 




RH 




HF 


7 


RHF 




. MC- 


1 


1 1 1 



4 



Matsubara Frequency \co^ 

Freericks, et al., Phys. Rev. B, Fig. 6 (c) 



24 



1.4 



1.2 - 



0.8 - 



0.6 



0.4 - 



0.2 - 







• 1 1 1 1 


1 1 1 1 

(a) 
H - 




- - RH - 




U IT - 

Mr 


-n 

\ 

•\\ 


RHF 

r\ n r 


" \\* 


• IVI ^ 




^ 4 • 1 • 4 m 



01 23456789 10 

Matsubara Frequency \co^ 

Freericks, et al., Phys. Rev. B, Fig. 7 (a) 



25 




Matsubara Frequency \(jJ^ 

Freericks, et aL, Phys. Rev. B, Fig. 7 (b) 



26 







-2 



-3 



1 ^ 


J— -• h» "I — -•— 

(c) 




H 


11 

It 


RH _ 


1 






HF 




RHF 




• MC 


• 

1 


1 1 1 



1 2 3 4 5 



Matsubara Frequency Icj^^ 

Freericks, et al., Phys. Rev. B, Fig. 7 (c) 



27 



0.4 







T — I — I — I — I — I — I — I — r 



H 



(a) 



- - RH 




01 23456789 10 



Matsubara Frequency \co 

Freericks, et al., Phys. Rev. B, Fig. 8 (a) 



28 



1 


(b) 




H 




RH 


— 1 


HF ~ 

1 1 1 


— 1 


RHF 




. MC 


1 


1 1 1 



1 2 3 4 5 

Matsubara Frequency \co^ 

Freericks, et al., Phys. Rev. B, Fig. 8 (b) 



29 



0.4 



0.3 - 



0.2 - 



0.1 - 









1 ll ' 

'/ ' 


1 1 






H 




l/ 


RH 






HF 




:// 






.7 


RHF 

1 \ 1 II 


— 


• MC - 
















f • 












L • 

1 r 1 


1 1 



0.2 0.4 0.6 0.8 1 

Interaction strength g/(t*+g) 

Freericks, et al., Phys. Rev. B, Fiq. 9 



30 




1 0.8 0.6 0.4 0.2 

Electron concentration 

Freericks, et aL, Phys. Rev. B, Fiq. 10 



31 



